この資料は「とりあえずTidyverseで出来ることを知る」という目的のもと、「使いそうだな」というTidyverseの関数をピックアップして紹介し、いくつかの実践的な例で処理を紹介する形式となっています。そのためTidyverseの設計思想や一連のパッケージの使い方を学ぶには不足があります。それらについてはR for Data Scienceのページで公式に詳しく解説されているので、そちらを見ていただければと思います。
またRそのものの使い方についてもかなり説明を省いている部分がありますので、それらについてはR-Tipsや関連書籍等で確認していただければと思います。
作図に関しても、ここでは「ggplot2を使ってとりあえず作図を行って、ある程度の体裁を整える」方法の解説を行っています。そのため、例としてして作成している図の体裁には、発表等に用いるにあたって不十分な部分があります。実際に作図を行われる際には雑誌の規定等を参照しながら各自で微調整を行ってください。
この資料で紹介する内容に基づいて起こる研究上の損害等については責任を負いかねますので、ここにあるコード等をご自身で用いられる場合には、公式のレファレンスでしっかり関数の仕様を理解する、自分自身でコードの確認を行うなどをした上で処理を実行してください。
R is a language and environment for statistical computing and graphics. It is a GNU project which is similar to the S language and environment which was developed at Bell Laboratories (formerly AT&T, now Lucent Technologies) by John Chambers and colleagues.
特徴:
※Pythonとの違いは? →個人的には、「Pythonはプログラミング言語、Rはコマンドで操作できるソフト」という印象を強く受ける。Rにおけるスクリプトは多くの場合、「プログラムを記述したもの」というよりも「データ処理の過程を記録したもの」である。
Rはコマンドで動きます。
基本的には、Rでの作業はオブジェクトとして格納されているデータを 関数に渡して、何かしらの処理をし、出力を得ることの繰り返しです。
この「データ」には、CSVファイルなどから読み込んだ外部データやR上で生成したデータ(乱数など)などがあります。
R_flow
R_merit
model <- glm(data = data, ...)等のコードからはそれらの情報を含む、解析方法についての十分な情報を得ることができるThe tidyverse is an opinionated collection of R packages designed for data science. All packages share an underlying design philosophy, grammar, and data structures.
Tidyverseはデータを扱うために一貫性をもって作成されたパッケージ群
Tidyverseではデータフレームに対する操作を関数で行うことができる
従来はデータフレームを入出力にできずに不便だった点も、Tidyverseでは解決されている
Ex. iris内における各種の観測数をカウント
table(iris$Species)
##
## setosa versicolor virginica
## 50 50 50
iris %>% count(Species)
## Species n
## 1 setosa 50
## 2 versicolor 50
## 3 virginica 50
素のRを使うよりも圧倒的に楽&わかりやすい(処理が明示的)
Ex.列を抜き出して、条件で選択
素のR
iris2 <- iris
iris2 <- iris2[, c("Sepal.Length", "Species")]
iris2 <- iris2[iris2[, "Sepal.Length"] < 5.0, ]
Tidyverse
iris2 <- iris %>%
select(Sepal.Length, Species) %>%
filter(Sepal.Length < 5.0)
生態学分野では多くの場合、調査で得たデータをデータフレームとしてRに読み込み、データ整形・解析・作図を行う。そのためデータフレーム志向で設計されたTidyverseは非常に使いやすく有用である(※個人的な意見です)
(簡単に言うと)処理を横につなげられるようにするツール
function(x, …) ⇔ x %>% function(…)
本当は段階的に別名オブジェクトに格納した方がコード的には確実だけど、 実際には連続で関数を適用したりするのでこっちの方が便利
unique(pull(filter(iris, Sepal.Length < 5.0), Sepal.Length))
## [1] 4.9 4.7 4.6 4.4 4.8 4.3 4.5
iris %>% filter(Sepal.Length < 5.0) %>% pull(Sepal.Length) %>% unique()
## [1] 4.9 4.7 4.6 4.4 4.8 4.3 4.5
パイプによってデータ処理のフローを、視覚的にもわかりやすく記述することができる
RStudio
RStudio is an integrated development environment (IDE) for R. It includes a console, syntax-highlighting editor that supports direct code execution, as well as tools for plotting, history, debugging and workspace management.
https://www.rstudio.com/products/rstudio/#rstudio-desktop より
RStudioはR用の統合解析環境である。
データや出力結果を置いておくためのディレクトリ(フォルダ)を指定する。
Rstudio上では Session > Set Working Directory > Choose Directory で設定できる。
今回はこちらのページでCode>Download ZIPと進み、ダウンロードして解凍したフォルダを作業ディレクトリにすることで以下コードを円滑に実行できる。
パッケージの読み込みはlibrary関数で行う。もしTidyverseがインストールされていない場合には、install.package("tidyverse")でパッケージをインストールする。
library(tidyverse)
# インストールは
# install.package("tidyverse")
今回は以下のデータペーパーで提供されているデータを例として使用する。
Schrader, J., Moeljono, S., Tambing, J., Sattler, C., & Kreft, H. 2020. A new dataset on plant occurrences on small islands, including species abundances and functional traits across different spatial scales. Biodiversity Data Journal 8:.
URL: https://bdj.pensoft.net/article/55275/instance/5664007/
Use license: Creative Commons Public Domain Waiver (CC-Zero)
インドネシアのRaja Ampat諸島に位置する60の小規模な島における木本植物のデータセット。 各島に2x2mのプロットからなるトランセクトを島の大きさに応じて設置し、プロット内に根を張る胸高直径2㎝以上のすべての種を記録したデータである。 その他に、島についての詳細データや植物種の機能的特性のデータ、種名のリスト、GBIF上に記録されている種の出現記録のデータが付属している。
oo_412824
Schrader et al. (2020) より
各データは以下のファイルに格納されている(元データを微調整し、CSV形式で保存したもの):
調査、データの詳細については詳しくはリンク先を参照のこと。
以下の項目名に示されている関数名にはレファレンスページへのリンクを貼ってあります。それぞれの関数を使う際には、ぜひ一度見てみてください(Rから呼び出せる関数ヘルプと同じ内容です)。
read_scv関数はread.csv関数とほぼ同等の関数であり、目的とするCSVファイル(やその他の形式のテキストファイル)のパスを渡すとデータの読み込みを行ってくれる。ただし、この際に読み込まれたデータはTibble形式で格納されることに注意。
※data.frameとtibbleは厳密には違いますが、この資料ではわかりやすさのために両方とも「データフレーム」と呼んで区別していません。
community <- read_csv("data/Community_data.csv")
island <- read_csv("data/Island_Data.csv")
trait <- read_csv("data/Plant_functional_trait_data.csv")
occurence <- read_csv("data/Plant_occurrences.csv")
species <- read_csv("data/Species_data.csv")
読み込んだデータの中身は以下のようになっている:
community
## # A tibble: 2,215 x 6
## island_ID transect_ID plot_ID species_ID DBH_cm tree_height_m
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 GB1 GB1_T2 GB1_T2_ST5 species_12 2 3
## 2 GB1 GB1_T5 GB1_T5_ST2 species_25 2 3.1
## 3 GB3 GB3_T1 GB3_T1_ST3 species_22 2 2.8
## 4 GB3 GB3_T3 GB3_T3_ST3 species_3 2 2.5
## 5 GB3 GB3_T3 GB3_T3_ST5 species_8 2 2.2
## 6 GB3 GB3_T4 GB3_T4_ST5 species_5 2 3.1
## 7 GB5 GB5_T1 GB5_T1_ST1 species_35 2 2.2
## 8 GB6 GB6_T2 GB6_T2_ST5 species_38 2 2.4
## 9 GB7 GB7_T1 GB7_T1_ST1 species_8 2 3
## 10 GB7 GB7_T1 GB7_T1_ST2 species_8 2 2.6
## # ... with 2,205 more rows
colnames(community)
## [1] "island_ID" "transect_ID" "plot_ID" "species_ID"
## [5] "DBH_cm" "tree_height_m"
island
## # A tibble: 60 x 12
## island_ID island_coordina~ island_area island_perimeter distance_Gam
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 GB1 -0.520461, 130.~ 4774. 303. 59.3
## 2 GB2 -0.517484, 130.~ 7.29 10.2 56.6
## 3 GB3 -0.517587, 130.~ 2330. 184. 172.
## 4 GB4 -0.517875, 130.~ 8.06 11.3 192.
## 5 GB5 -0.517319, 130.~ 20.3 16.6 136.
## 6 GB6 -0.515251, 130.~ 317. 85.1 382.
## 7 GB7 -0.515023, 130.~ 1575. 150. 345.
## 8 GB8 -0.516626, 130.~ 1264. 175. 107.
## 9 GB9 -0.514338, 130.~ 1716. 152. 400.
## 10 GB10 -0.516784, 130.~ 121. 46.4 115.
## # ... with 50 more rows, and 7 more variables: buffer_area_1000m <dbl>,
## # tree_basal_area <dbl>, species_number <dbl>, soil_depth_mean <dbl>,
## # leaf_litter_cover <dbl>, no_transects <dbl>, no_plots <dbl>
colnames(island)
## [1] "island_ID" "island_coordinates" "island_area"
## [4] "island_perimeter" "distance_Gam" "buffer_area_1000m"
## [7] "tree_basal_area" "species_number" "soil_depth_mean"
## [10] "leaf_litter_cover" "no_transects" "no_plots"
trait
## # A tibble: 57 x 13
## species_ID chlorophyll_SPAD chlorophyll_mod fruit_mass seed_mass LMA
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 species_1 50 59.2 0.07 14.4 0.91
## 2 species_2 47.4 54.7 0.22 30.4 1.94
## 3 species_3 59.2 77.3 0.14 41.5 1.23
## 4 species_4 46.0 52.5 0.13 0.51 1.47
## 5 species_5 54.4 67.6 NA 14.7 1.17
## 6 species_6 43.4 48.2 0.51 2.51 1.43
## 7 species_7 41.4 45.1 0.01 6 0.56
## 8 species_8 57.9 74.6 0.59 88.6 1.72
## 9 species_9 56.8 72.3 0.2 74.5 0.79
## 10 species_10 NA NA NA NA 1.17
## # ... with 47 more rows, and 7 more variables: leaf_area <dbl>,
## # wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## # leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>
colnames(trait)
## [1] "species_ID" "chlorophyll_SPAD" "chlorophyll_mod"
## [4] "fruit_mass" "seed_mass" "LMA"
## [7] "leaf_area" "wood_density" "max_tree_height"
## [10] "max_tree_height_3" "leaf_N" "leaf_C"
## [13] "leaf_P"
occurence
## # A tibble: 373 x 17
## id basisOfRecord occurrenceID recordedBy eventDate islandGroup country
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 2 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 3 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 4 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 5 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 6 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 7 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 8 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 9 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## 10 plan~ LivingSpecim~ plant-occur~ Julian Sc~ 2016-06-~ Raja Ampat Indone~
## # ... with 363 more rows, and 10 more variables: countryCode <chr>,
## # decimalLatitude <dbl>, decimalLongitude <dbl>, geodeticDatum <chr>,
## # coordinateUncertaintyInMeters <dbl>, identificationQualifier <chr>,
## # scientificName <chr>, kingdom <chr>, family <chr>, taxonRank <chr>
colnames(occurence)
## [1] "id" "basisOfRecord"
## [3] "occurrenceID" "recordedBy"
## [5] "eventDate" "islandGroup"
## [7] "country" "countryCode"
## [9] "decimalLatitude" "decimalLongitude"
## [11] "geodeticDatum" "coordinateUncertaintyInMeters"
## [13] "identificationQualifier" "scientificName"
## [15] "kingdom" "family"
## [17] "taxonRank"
species
## # A tibble: 57 x 5
## species_ID Family Species Author UNIPA_Voucher_ID
## <chr> <chr> <chr> <chr> <chr>
## 1 species_1 Rubiaceae Ixora timorens~ Decne. Schrader 12/13/82/99/~
## 2 species_2 Ebenaceae Diospyros mari~ Blume Schrader 14/30/40/80/~
## 3 species_3 Moraceae Ficus microcar~ L.f. Schrader 88/92/98/48/~
## 4 species_4 Moraceae Ficus tinctoria G.Forst. Schrader 11/97/155
## 5 species_5 Primulace~ Myrsine rawace~ A. DC. Schrader 35/95/173/101
## 6 species_6 Rubiaceae Timonius sp. 1 <NA> Schrader 20/151
## 7 species_7 Euphorbia~ Macaranga dioi~ (G.Forst.) M・l.~ Schrader 4/176/112
## 8 species_8 Myrtaceae Eugenia reinwa~ (Blume) DC. Schrader 10/121/159
## 9 species_9 Micromelum Micromelum min~ (G.Forst.) Wigh~ Schrader 106/147
## 10 species_10 Celastrac~ Pleurostylia o~ (Wall.) Alston Schrader 58
## # ... with 47 more rows
colnames(species)
## [1] "species_ID" "Family" "Species" "Author"
## [5] "UNIPA_Voucher_ID"
Tidy_data
整然データについて詳しくはこちらのページを参照のこと。
selectはデータフレームから特定の列を抜き出すために用いる関数である。この際、抜き出された列がたとえ一列だったとしてもデータフレームとして出力される。ベクトルとして抜き出したいのであれば、後述するpullを用いる
community %>%
select(island_ID, transect_ID, plot_ID, species_ID)
## # A tibble: 2,215 x 4
## island_ID transect_ID plot_ID species_ID
## <chr> <chr> <chr> <chr>
## 1 GB1 GB1_T2 GB1_T2_ST5 species_12
## 2 GB1 GB1_T5 GB1_T5_ST2 species_25
## 3 GB3 GB3_T1 GB3_T1_ST3 species_22
## 4 GB3 GB3_T3 GB3_T3_ST3 species_3
## 5 GB3 GB3_T3 GB3_T3_ST5 species_8
## 6 GB3 GB3_T4 GB3_T4_ST5 species_5
## 7 GB5 GB5_T1 GB5_T1_ST1 species_35
## 8 GB6 GB6_T2 GB6_T2_ST5 species_38
## 9 GB7 GB7_T1 GB7_T1_ST1 species_8
## 10 GB7 GB7_T1 GB7_T1_ST2 species_8
## # ... with 2,205 more rows
community %>%
select(-DBH_cm, -tree_height_m)
## # A tibble: 2,215 x 4
## island_ID transect_ID plot_ID species_ID
## <chr> <chr> <chr> <chr>
## 1 GB1 GB1_T2 GB1_T2_ST5 species_12
## 2 GB1 GB1_T5 GB1_T5_ST2 species_25
## 3 GB3 GB3_T1 GB3_T1_ST3 species_22
## 4 GB3 GB3_T3 GB3_T3_ST3 species_3
## 5 GB3 GB3_T3 GB3_T3_ST5 species_8
## 6 GB3 GB3_T4 GB3_T4_ST5 species_5
## 7 GB5 GB5_T1 GB5_T1_ST1 species_35
## 8 GB6 GB6_T2 GB6_T2_ST5 species_38
## 9 GB7 GB7_T1 GB7_T1_ST1 species_8
## 10 GB7 GB7_T1 GB7_T1_ST2 species_8
## # ... with 2,205 more rows
各列の値に基づく条件を用いて行を選択したい場合にはfilter関数を用いる。素のRで行う場合にはこの処理の記述は非常に煩雑になるため、filter関数は積極的に使っていきたい。
community %>%
filter((island_ID == "GB7")&(transect_ID == "GB7_T1")&(plot_ID == "GB7_T1_ST2"))
## # A tibble: 4 x 6
## island_ID transect_ID plot_ID species_ID DBH_cm tree_height_m
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 GB7 GB7_T1 GB7_T1_ST2 species_8 2 2.6
## 2 GB7 GB7_T1 GB7_T1_ST2 species_8 2.3 2.5
## 3 GB7 GB7_T1 GB7_T1_ST2 species_1 2.3 2
## 4 GB7 GB7_T1 GB7_T1_ST2 species_39 6.1 5.2
データフレーム中の列をベクトルとして抜き出したい場合にはpullを用いる。 データフレームに、新たなデータを含む新たな列を追加したい場合にはmutateを用いる。多くの場合、既存の列の値をもとに新たな値を算出し、その値を格納する処理に用いられる。 str_関数はパッケージstringrに含まれる関数であり、ここで用いているstr_extractは与えられた文字列から特定パターン(正規表現も可能)と一致する部分を抜き出してくれる関数である。
occurence %>% pull(scientificName) -> occurence_scientificName
occurence_scientificName[1:10]
## [1] "Glochidion castaneum Airy Shaw"
## [2] "Ixora timorensis Decne."
## [3] "Ficus pedunculosa Miq."
## [4] "Ficus prasinicarpa Elmer"
## [5] "Ficus microcarpa L.f."
## [6] "Ficus tinctoria G.Forst."
## [7] "Rapanea rawacensis (A. DC.) Mez"
## [8] "Aglaia elaeagnoidea (A.Juss.) Benth."
## [9] "Severinia lauterbachii Swingle"
## [10] "Planchonella obovata (R.Br.) Pierre"
occurence2 <- occurence %>%
mutate(species_name = str_extract(scientificName, pattern = "^[A-Za-z]+ [A-Za-z]+"))
occurence2 %>% pull(species_name) -> occurence2_scientificName
occurence2_scientificName[1:10]
## [1] "Glochidion castaneum" "Ixora timorensis" "Ficus pedunculosa"
## [4] "Ficus prasinicarpa" "Ficus microcarpa" "Ficus tinctoria"
## [7] "Rapanea rawacensis" "Aglaia elaeagnoidea" "Severinia lauterbachii"
## [10] "Planchonella obovata"
separateは特定列に含まれる文字列を、特定文字列を基準に分割し、2列に分けて格納する。ymdは日付を扱うためのパッケージlubridateに含まれる関数であり、「Year-Month-Day」形式で文字列として格納されている日付をDate型に変換してくれる。ここで用いている[パッケージ名]::[関数]という記法は特定のパッケージから1つだけ関数を呼び出したいときに便利。
occurence3 <- occurence2 %>%
separate(eventDate, sep = "/", into = c("Date1", "Date2")) %>%
mutate(
Date1 = lubridate::ymd(Date1),
Date2 = lubridate::ymd(Date2)
)
count関数は与えられた列内に含まれる各要素の観測数(=データフレーム内での行数)を数え、データフレーム形式で出力してくれる。2つ以上の列名を与えた場合には、各要素の組み合わせそれぞれの観測数を数えてくれる。 同様の処理はgroup_byとtallyを組み合わせることでも行うことができる。
community %>%
count(island_ID, species_ID)
## # A tibble: 390 x 3
## island_ID species_ID n
## <chr> <chr> <int>
## 1 GB1 species_12 1
## 2 GB1 species_2 3
## 3 GB1 species_22 2
## 4 GB1 species_25 7
## 5 GB1 species_26 2
## 6 GB1 species_32 1
## 7 GB1 species_37 5
## 8 GB1 species_38 5
## 9 GB1 species_40 3
## 10 GB1 species_43 1
## # ... with 380 more rows
group_byは与えられた列の要素に基づいてデータフレーム内の行をグループ化する関数である。summariseは各グループに含まれるデータに対し関数を適用することで値を要約し、結果をデータフレームとして返す。
community %>%
group_by(species_ID) %>%
summarise(
mean_DBH = mean(DBH_cm),
mean_height = mean(tree_height_m)
)
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 57 x 3
## species_ID mean_DBH mean_height
## <chr> <dbl> <dbl>
## 1 species_1 4.31 3.95
## 2 species_10 4.4 4.2
## 3 species_11 2.86 2.88
## 4 species_12 2.86 3.28
## 5 species_13 2.69 3.46
## 6 species_15 4.9 4.52
## 7 species_16 3.7 3.6
## 8 species_17 2.36 3.28
## 9 species_18 2.99 2.44
## 10 species_19 2.6 2
## # ... with 47 more rows
full_joinはbyで与えられた列をキーとして2つのデータフレームを結合する関数である。full_joinは2つのデータフレームのキー列に含まれるすべての要素を含む(つまりデータフレームAとBをfull_joinで結合した場合、データフレームAに記録がない要素についての行ではデータフレームAに含まれていた列の値がすべてNAとなる)。
このように、full_joinによる結合ではデータフレームの行数が気付かぬうちに増えてしまう危険性がある。そのため、適宜左側結合のleft_join、右側結合のright_join等と使い分けるべきである。
species
## # A tibble: 57 x 5
## species_ID Family Species Author UNIPA_Voucher_ID
## <chr> <chr> <chr> <chr> <chr>
## 1 species_1 Rubiaceae Ixora timorens~ Decne. Schrader 12/13/82/99/~
## 2 species_2 Ebenaceae Diospyros mari~ Blume Schrader 14/30/40/80/~
## 3 species_3 Moraceae Ficus microcar~ L.f. Schrader 88/92/98/48/~
## 4 species_4 Moraceae Ficus tinctoria G.Forst. Schrader 11/97/155
## 5 species_5 Primulace~ Myrsine rawace~ A. DC. Schrader 35/95/173/101
## 6 species_6 Rubiaceae Timonius sp. 1 <NA> Schrader 20/151
## 7 species_7 Euphorbia~ Macaranga dioi~ (G.Forst.) M・l.~ Schrader 4/176/112
## 8 species_8 Myrtaceae Eugenia reinwa~ (Blume) DC. Schrader 10/121/159
## 9 species_9 Micromelum Micromelum min~ (G.Forst.) Wigh~ Schrader 106/147
## 10 species_10 Celastrac~ Pleurostylia o~ (Wall.) Alston Schrader 58
## # ... with 47 more rows
trait
## # A tibble: 57 x 13
## species_ID chlorophyll_SPAD chlorophyll_mod fruit_mass seed_mass LMA
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 species_1 50 59.2 0.07 14.4 0.91
## 2 species_2 47.4 54.7 0.22 30.4 1.94
## 3 species_3 59.2 77.3 0.14 41.5 1.23
## 4 species_4 46.0 52.5 0.13 0.51 1.47
## 5 species_5 54.4 67.6 NA 14.7 1.17
## 6 species_6 43.4 48.2 0.51 2.51 1.43
## 7 species_7 41.4 45.1 0.01 6 0.56
## 8 species_8 57.9 74.6 0.59 88.6 1.72
## 9 species_9 56.8 72.3 0.2 74.5 0.79
## 10 species_10 NA NA NA NA 1.17
## # ... with 47 more rows, and 7 more variables: leaf_area <dbl>,
## # wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## # leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>
full_join(species, trait, by = c("species_ID" = "species_ID"))
## # A tibble: 57 x 17
## species_ID Family Species Author UNIPA_Voucher_ID chlorophyll_SPAD
## <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 species_1 Rubia~ Ixora ~ Decne. Schrader 12/13/~ 50
## 2 species_2 Ebena~ Diospy~ Blume Schrader 14/30/~ 47.4
## 3 species_3 Morac~ Ficus ~ L.f. Schrader 88/92/~ 59.2
## 4 species_4 Morac~ Ficus ~ G.For~ Schrader 11/97/~ 46.0
## 5 species_5 Primu~ Myrsin~ A. DC. Schrader 35/95/~ 54.4
## 6 species_6 Rubia~ Timoni~ <NA> Schrader 20/151 43.4
## 7 species_7 Eupho~ Macara~ (G.Fo~ Schrader 4/176/~ 41.4
## 8 species_8 Myrta~ Eugeni~ (Blum~ Schrader 10/121~ 57.9
## 9 species_9 Micro~ Microm~ (G.Fo~ Schrader 106/147 56.8
## 10 species_10 Celas~ Pleuro~ (Wall~ Schrader 58 NA
## # ... with 47 more rows, and 11 more variables: chlorophyll_mod <dbl>,
## # fruit_mass <dbl>, seed_mass <dbl>, LMA <dbl>, leaf_area <dbl>,
## # wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## # leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>
pivot_longerはデータフレームを縦長変換するための関数であり、与えられた列の列名を一つの列に、各列の値をもう一つの列に格納する。特に読み込んだデータフレームが非Tidyな場合(例えば、各列に各月の値が入っている場合など)に、それらをTidyなデータに変換するために用いることができる。
island_pivot1 <- island %>%
select(island_ID, island_coordinates, island_area, island_perimeter)
island_pivot2 <- island_pivot1 %>%
pivot_longer(cols = c(island_area, island_perimeter), names_to = "Island_property", values_to = "Value")
island_pivot2
## # A tibble: 120 x 4
## island_ID island_coordinates Island_property Value
## <chr> <chr> <chr> <dbl>
## 1 GB1 -0.520461, 130.580800 island_area 4774.
## 2 GB1 -0.520461, 130.580800 island_perimeter 303.
## 3 GB2 -0.517484, 130.575020 island_area 7.29
## 4 GB2 -0.517484, 130.575020 island_perimeter 10.2
## 5 GB3 -0.517587, 130.568534 island_area 2330.
## 6 GB3 -0.517587, 130.568534 island_perimeter 184.
## 7 GB4 -0.517875, 130.568515 island_area 8.06
## 8 GB4 -0.517875, 130.568515 island_perimeter 11.3
## 9 GB5 -0.517319, 130.569511 island_area 20.3
## 10 GB5 -0.517319, 130.569511 island_perimeter 16.6
## # ... with 110 more rows
pivot_widerはデータフレームの横長変換を行うための関数であり、一つの列に含まれる要素を列名、もう一つの列に含まれる値を各列の値として引き出しデータフレームを横長に展開する。横長変換する際にはデータが存在しないセルも生じるが、その際にはそのセルの値はNAとなる。
pivot_longerの逆にあたる処理を行う関数。
island_pivot3 <- island_pivot2 %>%
pivot_wider(names_from = "Island_property", values_from = "Value")
island_pivot3
## # A tibble: 60 x 4
## island_ID island_coordinates island_area island_perimeter
## <chr> <chr> <dbl> <dbl>
## 1 GB1 -0.520461, 130.580800 4774. 303.
## 2 GB2 -0.517484, 130.575020 7.29 10.2
## 3 GB3 -0.517587, 130.568534 2330. 184.
## 4 GB4 -0.517875, 130.568515 8.06 11.3
## 5 GB5 -0.517319, 130.569511 20.3 16.6
## 6 GB6 -0.515251, 130.568484 317. 85.1
## 7 GB7 -0.515023, 130.569425 1575. 150.
## 8 GB8 -0.516626, 130.572132 1264. 175.
## 9 GB9 -0.514338, 130.570002 1716. 152.
## 10 GB10 -0.516784, 130.571887 121. 46.4
## # ... with 50 more rows
write_csvはwrite.csvとほぼ同等の関数で、データフレームを指定されたファイル名のテキストファイル(CSV含む)に出力する。デフォルトのエンコーディングはUTF-8であるため、日本語が含まれる場合は注意。
write_csv(occurence3, file = "occurence3.csv")
ここまで見て来たTidyverseの関数を使って以下のような実践的な処理をどのように行うことができるか、考えてみよう。
紹介していない関数も使うので、適宜検索しながら進めよう。
種ごとのDBHの平均と標準誤差、種ごとの樹高の平均と標準誤差と種のデータ(species、trait)を結合したデータフレームを作成する
島ごとの総出現種数、総出現属数、総出現科数と島のデータを結合したデータフレームを作成する
調査プロット名×種名行列(=PCAやDCAで使える形)を作成する
答えを見る前に考えてみてね
group_byとsummariseを用いて、各種ごとに平均値に加えて標準誤差を算出する。得られたデータフレームをspecies_IDをキーとしてspecies、traitとfull_joinする。
community %>%
group_by(species_ID) %>%
summarise(
mean_DBH = mean(DBH_cm),
se_DBH = sd(DBH_cm)/sqrt(length(DBH_cm)),
mean_height = mean(tree_height_m),
se_height = sd(tree_height_m)/sqrt(length(tree_height_m))
) -> species_dbh_height
## `summarise()` ungrouping output (override with `.groups` argument)
species_perperty <- full_join(species, trait, by = "species_ID") %>% full_join(species_dbh_height, by = "species_ID")
species_perperty
## # A tibble: 57 x 21
## species_ID Family Species Author UNIPA_Voucher_ID chlorophyll_SPAD
## <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 species_1 Rubia~ Ixora ~ Decne. Schrader 12/13/~ 50
## 2 species_2 Ebena~ Diospy~ Blume Schrader 14/30/~ 47.4
## 3 species_3 Morac~ Ficus ~ L.f. Schrader 88/92/~ 59.2
## 4 species_4 Morac~ Ficus ~ G.For~ Schrader 11/97/~ 46.0
## 5 species_5 Primu~ Myrsin~ A. DC. Schrader 35/95/~ 54.4
## 6 species_6 Rubia~ Timoni~ <NA> Schrader 20/151 43.4
## 7 species_7 Eupho~ Macara~ (G.Fo~ Schrader 4/176/~ 41.4
## 8 species_8 Myrta~ Eugeni~ (Blum~ Schrader 10/121~ 57.9
## 9 species_9 Micro~ Microm~ (G.Fo~ Schrader 106/147 56.8
## 10 species_10 Celas~ Pleuro~ (Wall~ Schrader 58 NA
## # ... with 47 more rows, and 15 more variables: chlorophyll_mod <dbl>,
## # fruit_mass <dbl>, seed_mass <dbl>, LMA <dbl>, leaf_area <dbl>,
## # wood_density <dbl>, max_tree_height <dbl>, max_tree_height_3 <dbl>,
## # leaf_N <dbl>, leaf_C <dbl>, leaf_P <dbl>, mean_DBH <dbl>, se_DBH <dbl>,
## # mean_height <dbl>, se_height <dbl>
以下のようにselect、distinct、countを組み合わせることで、各島に存在する木本植物の種数、属数、科数を算出し、最後にisland_IDをキーに結合する。 属数に関しては、str_extractで種名から属名を抽出し、それを参照してカウントする。
ただし、60島のデータを含むislandとそのうちの40島における調査データをもとに作成したデータフレームをfull_joinで結合しているため、island_propertyにはn_species等の列の値がNAになっている20行のデータが含まれることになる点には注意。
community_species <- left_join(community, species, by = "species_ID")
sp_per_island <- community_species %>%
select(island_ID, species_ID) %>%
distinct() %>%
count(island_ID, name = "n_species")
sp_per_island
## # A tibble: 40 x 2
## island_ID n_species
## <chr> <int>
## 1 GB1 18
## 2 GB10 5
## 3 GB11 9
## 4 GB12 14
## 5 GB13 10
## 6 GB14 11
## 7 GB15 9
## 8 GB16 6
## 9 GB17 3
## 10 GB18 8
## # ... with 30 more rows
genera_per_island <- community_species %>%
mutate(Genera = str_extract(Species, pattern = "[A-Za-z]+")) %>%
select(island_ID, Genera) %>%
distinct() %>%
count(island_ID, name = "n_genera")
genera_per_island
## # A tibble: 40 x 2
## island_ID n_genera
## <chr> <int>
## 1 GB1 17
## 2 GB10 5
## 3 GB11 9
## 4 GB12 13
## 5 GB13 10
## 6 GB14 11
## 7 GB15 9
## 8 GB16 6
## 9 GB17 3
## 10 GB18 7
## # ... with 30 more rows
family_per_island <- community_species %>%
count(island_ID, Family) %>%
select(island_ID, Family) %>%
distinct() %>%
count(island_ID, name = "n_family")
family_per_island
## # A tibble: 40 x 2
## island_ID n_family
## <chr> <int>
## 1 GB1 14
## 2 GB10 5
## 3 GB11 6
## 4 GB12 9
## 5 GB13 8
## 6 GB14 10
## 7 GB15 8
## 8 GB16 6
## 9 GB17 2
## 10 GB18 7
## # ... with 30 more rows
island_property <- full_join(island, sp_per_island, by = "island_ID") %>%
full_join(genera_per_island, by = "island_ID") %>%
full_join(family_per_island, by = "island_ID")
island_property
## # A tibble: 60 x 15
## island_ID island_coordina~ island_area island_perimeter distance_Gam
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 GB1 -0.520461, 130.~ 4774. 303. 59.3
## 2 GB2 -0.517484, 130.~ 7.29 10.2 56.6
## 3 GB3 -0.517587, 130.~ 2330. 184. 172.
## 4 GB4 -0.517875, 130.~ 8.06 11.3 192.
## 5 GB5 -0.517319, 130.~ 20.3 16.6 136.
## 6 GB6 -0.515251, 130.~ 317. 85.1 382.
## 7 GB7 -0.515023, 130.~ 1575. 150. 345.
## 8 GB8 -0.516626, 130.~ 1264. 175. 107.
## 9 GB9 -0.514338, 130.~ 1716. 152. 400.
## 10 GB10 -0.516784, 130.~ 121. 46.4 115.
## # ... with 50 more rows, and 10 more variables: buffer_area_1000m <dbl>,
## # tree_basal_area <dbl>, species_number <dbl>, soil_depth_mean <dbl>,
## # leaf_litter_cover <dbl>, no_transects <dbl>, no_plots <dbl>,
## # n_species <int>, n_genera <int>, n_family <int>
countで調査プロット×植物種のクロス集計を行い、得られたデータフレームをpivot_widerで横長返還してからas.matrixで行列に変換する。 その際、pivot_widerを行って生成されたNAはmutate_allとreplace_naで0に置換し、列として格納されていたplot_IDはcolumn_to_rownamesで行名に移動する。
community %>%
count(plot_ID, species_ID) %>%
pivot_wider(names_from = species_ID, values_from = n) %>%
mutate_all(replace_na, 0) %>%
column_to_rownames("plot_ID") %>%
as.matrix() -> species_matrix
species_matrix[1:10, 1:10]
## species_5 species_8 species_80 species_82 species_22 species_12
## GB1_T1_ST1 1 1 1 0 0 0
## GB1_T1_ST2 2 5 0 0 0 0
## GB1_T1_ST3 0 6 0 0 0 0
## GB1_T1_ST4 0 2 1 0 0 0
## GB1_T1_ST5 0 1 0 0 0 0
## GB1_T2_ST1 0 2 0 1 0 0
## GB1_T2_ST2 0 2 0 0 1 0
## GB1_T2_ST4 0 1 0 0 0 0
## GB1_T2_ST5 0 3 0 0 0 1
## GB1_T3_ST1 0 1 0 0 0 0
## species_2 species_63 species_26 species_43
## GB1_T1_ST1 0 0 0 0
## GB1_T1_ST2 0 0 0 0
## GB1_T1_ST3 0 0 0 0
## GB1_T1_ST4 0 0 0 0
## GB1_T1_ST5 0 0 0 0
## GB1_T2_ST1 0 0 0 0
## GB1_T2_ST2 0 0 0 0
## GB1_T2_ST4 0 0 0 0
## GB1_T2_ST5 0 0 0 0
## GB1_T3_ST1 1 4 0 0
楽&わかりやすい(レイヤーが明示的)&やり直しやすい
plot(iris$Sepal.Length, iris$Sepal.Width)
# PNGで保存
png("iris_scatterplot.png", width = 480, height = 480)
plot(iris$Sepal.Length, iris$Sepal.Width)
dev.off()
## png
## 2
g <- ggplot() +
geom_point(data = iris, aes(x = Sepal.Length, y = Sepal.Width, color = Species))
plot(g)
# PNGで保存
ggsave(g, file = "iris_scatterggplot.png", width = 7, height = 7)
ggplot_layer
ggplot2における作図はggplot関数で元レイヤを作り、geom_*関数でデータをマッピングするところから始まる。 使う機会が非常に多いであろうgeom_point関数は与えられたデータを元に点をプロットする。
先に述べたように、island_propertyの60行中20行ではn_speciesの値がNAとなっている。よって、描画時にはこれらの値は描画されず、Removed 20 rows containing missing values (geom_point).という警告が出ることになる。
ggplot(data = island_property) +
geom_point(aes(x = island_area, y = n_species))
## Warning: Removed 20 rows containing missing values (geom_point).
g_point <- ggplot(data = island_property) +
geom_point(aes(x = island_area, y = n_species))
ggsave(g_point, file = "g_point.png")
## Saving 7 x 5 in image
## Warning: Removed 20 rows containing missing values (geom_point).
geom_point1
geom_*関数による各プロットは、レイヤとして重ね合わせることができる。
ggplot(data = island_property) +
geom_point(aes(x = island_area, y = n_species), color = "steelblue")
## Warning: Removed 20 rows containing missing values (geom_point).
ggplot(data = island_property) +
geom_point(aes(x = island_area, y = n_genera), color = "lightcoral")
## Warning: Removed 20 rows containing missing values (geom_point).
ggplot(data = island_property) +
geom_point(aes(x = island_area, y = n_family), color = "forestgreen")
## Warning: Removed 20 rows containing missing values (geom_point).
ggplot(data = island_property) +
geom_point(aes(x = island_area, y = n_species), color = "steelblue") +
geom_point(aes(x = island_area, y = n_genera), color = "lightcoral") +
geom_point(aes(x = island_area, y = n_family), color = "forestgreen")
## Warning: Removed 20 rows containing missing values (geom_point).
## Warning: Removed 20 rows containing missing values (geom_point).
## Warning: Removed 20 rows containing missing values (geom_point).
geom_pointに渡すaes関数内でcolor = <列名>を指定することにより、その列に含まれる値に基づいてプロット上の点の色分けを行うことができる。列に含まれる値が文字列やfactor含む離散値だった場合には色「分け」され、連続値だった場合には数値に応じたグラデーションによる描き分けがなされる。
※ここでは例示のためisland_property2を描画しているが、縦に3つ並ぶ値はすべて同じ島についてのものであることに注意。本来であればきちんとした整然データを描画することが望ましい(のだけれど、良い感じのカテゴリ値がなかったのでご容赦を)。
island_property2 <- island_property %>%
pivot_longer(cols = c(n_species, n_genera, n_family), names_to = "hierarchy", values_to = "value")
ggplot(data = island_property2) +
geom_point(aes(x = island_area, y = value, color = hierarchy)) +
scale_color_manual(values = c("forestgreen", "lightcoral", "steelblue"))
## Warning: Removed 60 rows containing missing values (geom_point).
geom_point3
ヒストグラムの描画はgeom_histogramで行う。aesにはx軸として指定したい列を与える。階級幅は、binwidthやbreaks等のパラメータで設定する。
ggplot(data = community) +
geom_histogram(aes(x = tree_height_m), breaks = seq(0, 18, 1))
箱ひげ図の描画はgeom_boxplotで行う。
ggplot(data = community) +
geom_boxplot(aes(x = island_ID, y = tree_height_m))
棒グラフの描画にはgeom_barを用いる。y軸の値をデータフレーム内から参照する場合にはstat = "identity"を、x軸とした値それぞれの観測数にする場合にはstat = "count"を指定する。また、積み上げ棒グラフにしたい場合にはaes内でfill = <分けたい要素が格納された列名>を指定する。割合を示したい際には、position = "fill"を指定する。
ggplot(data = community_species) +
geom_bar(aes(x = island_ID, fill = Family), stat = "count")
scale_*関数を用いることで、x軸やy軸の目盛設定や軸ラベル設定、色分けの詳細設定などを行うことができる。
ggplot(data = island_property2) +
geom_point(aes(x = island_area, y = value, color = hierarchy, shape = hierarchy), alpha = 0.7) +
scale_x_continuous(breaks = seq(0, 12000, by = 2000)) +
scale_y_continuous(breaks = seq(0, 30, by = 10)) +
scale_color_manual(name = "Classification level", values = c("forestgreen", "lightcoral", "steelblue"), labels = c("Family", "Genera", "Species")) +
scale_shape_discrete(name = "Classification level", labels = c("Family", "Genera", "Species")) +
xlab("Island area") +
ylab("Number of taxa")
## Warning: Removed 60 rows containing missing values (geom_point).
other_settings
geom_関数で描画を行い、scale_関数で色分け等の微調整を行ったとしても、最終的にPNG形式等で出力した図の体裁が整うわけではない。 パネル領域の背景、軸ラベルの文字、凡例の位置等の「見た目上の調整」はtheme関数内で行う必要がある。 ここで注意しなければならないのが「文字や点の大きさ」であり、これらは出力画像サイズを決定してから調整することをおすすめする。
g_point_color2_theme <- ggplot(data = island_property2) +
geom_point(aes(x = island_area, y = value, color = hierarchy, shape = hierarchy), size = 6, alpha = 0.7) +
scale_x_continuous(breaks = seq(0, 12000, by = 2000)) +
scale_y_continuous(breaks = seq(0, 30, by = 10)) +
scale_color_manual(name = "Classification level", values = c("forestgreen", "lightcoral", "steelblue"), labels = c("Family", "Genera", "Species")) +
scale_shape_discrete(name = "Classification level", labels = c("Family", "Genera", "Species")) +
xlab("Island area") +
ylab("Number of taxa") +
theme(
text = element_text(family = "sans"),
axis.title = element_text(size = 30, color = "black"),
axis.text = element_text(size = 30, color = "black"),
axis.line = element_line(color = "black"),
axis.ticks.x = element_line(color = "black"),
axis.ticks.length.x = unit(0.3, "cm"),
axis.ticks.y = element_line(color = "black"),
axis.ticks.length.y = unit(0.3, "cm"),
panel.border = element_rect(fill = NA, color = "black"),
plot.margin = margin(2, 2, 2, 2, "cm"),
legend.background = element_blank(),
legend.position = c(0.80, 0.15),
legend.title = element_text(size = 30),
legend.text = element_text(size = 30),
legend.box.background = element_blank()
)
ggsave(g_point_color2_theme, file = "g_point_color2_theme.png", width = 14, height = 14, dpi = 200)
## Warning: Removed 60 rows containing missing values (geom_point).
g_point_color2_theme
theme1
theme設定は出力画像サイズがほぼ変わらないのであれば、そのまま他のプロットでも使い回すことができる。これにより、複数の図のレイアウトに統一性を持たせることができる。
g_hist_theme <- ggplot(data = community) +
geom_histogram(aes(x = tree_height_m), breaks = seq(0, 18, 1), fill = "forestgreen") +
theme(
text = element_text(family = "sans"),
axis.title = element_text(size = 30, color = "black"),
axis.text = element_text(size = 30, color = "black"),
axis.line = element_line(color = "black"),
axis.ticks.x = element_line(color = "black"),
axis.ticks.length.x = unit(0.3, "cm"),
axis.ticks.y = element_line(color = "black"),
axis.ticks.length.y = unit(0.3, "cm"),
panel.border = element_rect(fill = NA, color = "black"),
plot.margin = margin(2, 2, 2, 2, "cm"),
legend.background = element_blank(),
legend.position = c(0.80, 0.15),
legend.title = element_text(size = 30),
legend.text = element_text(size = 30),
legend.box.background = element_blank()
)
ggsave(g_hist_theme, file = "g_hist_theme.png", width = 14, height = 14, dpi = 200)
g_hist_theme
theme設定はオブジェクトとして保存して呼び出すことができる。またすでにthemeで定義されている要素に関しても、あとからthemeで上書き再設定を行うことができる。
theme_example <- theme(
text = element_text(family = "sans"),
axis.title = element_text(size = 30, color = "black"),
axis.text = element_text(size = 30, color = "black"),
axis.line = element_line(color = "black"),
axis.ticks.x = element_line(color = "black"),
axis.ticks.length.x = unit(0.3, "cm"),
axis.ticks.y = element_line(color = "black"),
axis.ticks.length.y = unit(0.3, "cm"),
panel.border = element_rect(fill = NA, color = "black"),
plot.margin = margin(2, 2, 2, 2, "cm"),
legend.background = element_blank(),
legend.position = c(0.80, 0.15),
legend.title = element_text(size = 30),
legend.text = element_text(size = 30),
legend.box.background = element_blank()
)
g_hist_theme_overwrite <- ggplot(data = community) +
geom_histogram(aes(x = tree_height_m), breaks = seq(0, 18, 1), fill = "forestgreen") +
theme_example +
theme(
panel.background = element_blank(),
panel.grid.major = element_line(color = "black", linetype = 2)
)
ggsave(g_hist_theme_overwrite, file = "g_hist_theme_overwrite.png", width = 14, height = 14, dpi = 200)
g_hist_theme_overwrite
g_community_species <- ggplot(data = community_species) +
geom_bar(aes(x = island_ID, fill = Family), stat = "count") +
xlab("Island ID") +
ylab("Number of observations") +
theme_example +
theme(
axis.text.x = element_text(angle = -45, hjust = 0),
legend.position = c(0.82, 0.70)
)
ggsave(g_community_species, file = "g_community_species.png", width = 24, height = 14, dpi = 200)
g_community_species
species_propertyに含まれるデータをもとに、各種の平均DBHと平均樹高をプロットし、科名で色分けする。ここではそれに加えて、geom_segmentを用いて標準誤差をもとに95%信頼区間も表示している。
g_dbh_and_height <- ggplot(data = species_perperty) +
geom_point(aes(x = mean_DBH, y = mean_height, color = Family), size = 4) +
geom_segment(aes(x = mean_DBH, xend = mean_DBH, y = mean_height - 1.96*se_height, yend = mean_height + 1.96*se_height, color = Family)) +
geom_segment(aes(x = mean_DBH - 1.96*se_DBH, xend = mean_DBH + 1.96*se_DBH, y = mean_height, yend = mean_height, color = Family)) +
scale_x_continuous(name = "Average DBH (cm)", limits = c(0, 25)) +
scale_y_continuous(name = "Average height (m)", limits = c(0, 15)) +
scale_color_discrete(name = "Family of speceis") +
coord_cartesian(expand = FALSE) +
theme_example +
theme(
legend.position = c(0.80, 0.20),
legend.text = element_text(size = 15)
)
ggsave(g_dbh_and_height, file = "g_dbh_and_height.png", width = 14, height = 14, dpi = 200)
## Warning: Removed 4 rows containing missing values (geom_segment).
## Warning: Removed 4 rows containing missing values (geom_segment).
g_dbh_and_height
大きなデータを解析していると複数の同形式のデータについて同様の処理を適用したいケースが出てくる
具体的には、地点ごとに生物多様性指数を算出したい、種ごとのデータに対して回帰分析を行いたい場合などが考えられる
「基礎編」で扱ったグループごとに平均値を算出するような場合はいいが、処理が複雑になった場合にも対処できる方法を知らないと、同じコードを何回もコピペしてそれぞれを訂正するという非常に危ないコーディングを行うことになる
Rの基本パッケージを用いる場合、リスト化したデータとlapplyを用いることなどが方法の一つである
iris_sp <- list(Set = iris[iris$Species == "setosa", ], ver = iris[iris$Species == "versicolor", ], vir = iris[iris$Species == "virginica", ])
correlation_sp <- lapply(iris_sp, function(x) cor.test(x$Sepal.Length, x$Sepal.Width, method = "pearson"))
correlation_sp[[1]]
##
## Pearson's product-moment correlation
##
## data: x$Sepal.Length and x$Sepal.Width
## t = 7.6807, df = 48, p-value = 6.71e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5851391 0.8460314
## sample estimates:
## cor
## 0.7425467
Tidyverseではtidyr::nestとpurrr::mapを用いることでより簡単に、明示的に同様の処理を行い、結果をそのままデータフレーム(Tibble)に格納することができる。
iris %>%
group_by(Species) %>%
nest() %>%
mutate(correlation_sp = map(data, function(x) cor.test(x$Sepal.Length, x$Sepal.Width, method = "pearson")))
## # A tibble: 3 x 3
## # Groups: Species [3]
## Species data correlation_sp
## <fct> <list> <list>
## 1 setosa <tibble [50 x 4]> <htest>
## 2 versicolor <tibble [50 x 4]> <htest>
## 3 virginica <tibble [50 x 4]> <htest>
nest_map
今回のデータで島ごとに各種の観測数を集計したい場合には次のような方法が考えられる。ここで島ごとのリストはdata2の各行の値として格納されている。
community %>%
group_by(island_ID) %>%
nest() %>%
mutate(
data2 = map(data, function(data) {
count_of_species <- data %>%
group_by(species_ID) %>%
tally()
return(count_of_species)
}
)
) -> sp_count_per_island
sp_count_per_island
## # A tibble: 40 x 3
## # Groups: island_ID [40]
## island_ID data data2
## <chr> <list> <list>
## 1 GB1 <tibble [92 x 5]> <tibble [18 x 2]>
## 2 GB3 <tibble [84 x 5]> <tibble [17 x 2]>
## 3 GB5 <tibble [3 x 5]> <tibble [2 x 2]>
## 4 GB6 <tibble [45 x 5]> <tibble [8 x 2]>
## 5 GB7 <tibble [70 x 5]> <tibble [12 x 2]>
## 6 GB8 <tibble [72 x 5]> <tibble [13 x 2]>
## 7 GB9 <tibble [91 x 5]> <tibble [18 x 2]>
## 8 GB10 <tibble [11 x 5]> <tibble [5 x 2]>
## 9 GB11 <tibble [77 x 5]> <tibble [9 x 2]>
## 10 GB12 <tibble [62 x 5]> <tibble [14 x 2]>
## # ... with 30 more rows
sp_count_per_island[9, 3] %>% unnest()
## Warning: `cols` is now required when using unnest().
## Please use `cols = c(data2)`
## # A tibble: 9 x 2
## species_ID n
## <chr> <int>
## 1 species_1 10
## 2 species_2 2
## 3 species_23 6
## 4 species_38 8
## 5 species_39 1
## 6 species_41 1
## 7 species_5 6
## 8 species_61 2
## 9 species_8 41
count_per_island
島ごとに各種の観測数を算出してしまえば、それをもとに多様性の指標であるShannon’s Hを算出することができる。 ここでは種iの観測数を島内における全種の観測数で除したものをpiとして便宜的に用いている。
※「H’」の表記もありうるが、これも便宜的に(コードを書く都合上)「H」としている。
sp_count_per_island %>%
mutate(
Shannon_H = map_dbl(data2, function(data2){
percentage <- data2$n / sum(data2$n)
H = - sum(sapply(percentage, function(x) {x * log(x)}))
return(H)
}
)
) -> diversity_per_island
diversity_per_island
## # A tibble: 40 x 4
## # Groups: island_ID [40]
## island_ID data data2 Shannon_H
## <chr> <list> <list> <dbl>
## 1 GB1 <tibble [92 x 5]> <tibble [18 x 2]> 2.40
## 2 GB3 <tibble [84 x 5]> <tibble [17 x 2]> 2.37
## 3 GB5 <tibble [3 x 5]> <tibble [2 x 2]> 0.637
## 4 GB6 <tibble [45 x 5]> <tibble [8 x 2]> 1.14
## 5 GB7 <tibble [70 x 5]> <tibble [12 x 2]> 1.91
## 6 GB8 <tibble [72 x 5]> <tibble [13 x 2]> 2.19
## 7 GB9 <tibble [91 x 5]> <tibble [18 x 2]> 2.25
## 8 GB10 <tibble [11 x 5]> <tibble [5 x 2]> 1.39
## 9 GB11 <tibble [77 x 5]> <tibble [9 x 2]> 1.54
## 10 GB12 <tibble [62 x 5]> <tibble [14 x 2]> 2.08
## # ... with 30 more rows
diversity_per_island
island_diversity <- full_join(sp_per_island, diversity_per_island, by = "island_ID") %>% full_join(island, by = "island_ID")
g_area_diversity <- ggplot(data = island_diversity) +
geom_point(aes(x = island_area, y = Shannon_H), size = 6, alpha = 0.5) +
scale_x_continuous(name = "Island area", breaks = seq(0, 12000, by = 2000)) +
scale_y_continuous(name = "Shannon's H", breaks = seq(0, 3, by = 1)) +
theme_example
ggsave(g_area_diversity, file = "g_area_diversity.png", width = 14, height = 14, dpi = 200)
## Warning: Removed 20 rows containing missing values (geom_point).
g_distance_diversity <- ggplot(data = island_diversity) +
geom_point(aes(x = distance_Gam, y = Shannon_H), size = 6, alpha = 0.5) +
scale_x_continuous(name = "Shortest distance of each island to the nearest large landmass", breaks = seq(0, 1400, by = 200)) +
scale_y_continuous(name = "Shannon's H", breaks = seq(0, 3, by = 1)) +
theme_example
ggsave(g_distance_diversity, file = "g_distance_diversity.png", width = 14, height = 14, dpi = 200)
## Warning: Removed 20 rows containing missing values (geom_point).
今回扱っているデータのように、調査地点の座標情報がわかっている場合にはデータ解析の結果を主題図のように地図上に描画することができる
read_sf関数を用いることで、.shpや.geojesonで保存されているGISデータをsfオブジェクトとして読み込むことができる。 既存のデータフレームの座標情報を参照してそのデータをsfオブジェクトとしたい場合には、st_as_sf関数を用いる。その際には座標系をcrsに指定する必要があることに注意。
今回はベースマップとしてこちらのページで提供されているOpenStreetMapの陸域データをクリップしたものを使用した。
## install.pacages(sf)
library(sf)
## Warning: package 'sf' was built under R version 4.0.5
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
basemap <- read_sf("data/basemap_osm.geojson")
island_diversity_sf <- island_diversity %>%
separate(island_coordinates, sep = ", ", into = c("LAT", "LON")) %>%
mutate(LAT = str_replace(LAT, pattern = ",", replacement = "."), LON = str_replace(LON, pattern = ",", replacement = ".")) %>%
mutate(LAT = as.numeric(LAT), LON = as.numeric(LON)) %>%
st_as_sf(coords = c("LON", "LAT"), crs = "wgs84")
island_diversity_sf
## Simple feature collection with 60 features and 15 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 130.5577 ymin: -0.522189 xmax: 130.594 ymax: -0.505713
## Geodetic CRS: WGS 84
## # A tibble: 60 x 16
## island_ID n_species data data2 Shannon_H island_area island_perimeter
## * <chr> <int> <lis> <lis> <dbl> <dbl> <dbl>
## 1 GB1 18 <tib~ <tib~ 2.40 4774. 303.
## 2 GB10 5 <tib~ <tib~ 1.39 121. 46.4
## 3 GB11 9 <tib~ <tib~ 1.54 817. 120.
## 4 GB12 14 <tib~ <tib~ 2.08 1650. 154.
## 5 GB13 10 <tib~ <tib~ 2.12 602. 90.0
## 6 GB14 11 <tib~ <tib~ 1.99 535. 90.2
## 7 GB15 9 <tib~ <tib~ 1.58 381. 77.8
## 8 GB16 6 <tib~ <tib~ 1.20 137. 43.5
## 9 GB17 3 <tib~ <tib~ 0.509 18.4 18.1
## 10 GB18 8 <tib~ <tib~ 1.78 433. 77
## # ... with 50 more rows, and 9 more variables: distance_Gam <dbl>,
## # buffer_area_1000m <dbl>, tree_basal_area <dbl>, species_number <dbl>,
## # soil_depth_mean <dbl>, leaf_litter_cover <dbl>, no_transects <dbl>,
## # no_plots <dbl>, geometry <POINT [arc_degree]>
sfオブジェクトはgeom_sf関数によってggplot2レイヤーとして描画することができる。この際、座標情報はsfオブジェクトから自動的に呼び出されるため、aes内では点やポリゴンの色分けに使う列名等を指定することになる。 また描画範囲はcoord_sf関数のxlimおよびylimで指定できる。
※今回の場合、island_propertyには多様性指数を算出できない(=調査記録がない)20島のデータが含まれているため、それらの点は灰色の点となっている。
g_gis_raw <- ggplot() +
geom_sf(data = basemap, color = "black", fill = "darkolivegreen3") +
geom_sf(data = island_diversity_sf, aes(color = Shannon_H), size = 4, alpha = 0.8) +
scale_color_viridis_c(name = "Shannon's H", breaks = seq(0, 3, by = 1)) +
coord_sf(xlim = c(130.55, 130.60), ylim = c(-0.53, -0.5)) +
labs(caption = "Land data from OSM(https://osmdata.openstreetmap.de/data/land-polygons.html)") +
theme_example +
theme(
legend.position = c(0.15, 0.20)
)
ggsave(g_gis_raw, file = "g_gis_raw.png", width = 14, height = 8, dpi = 200)
g_gis_raw
sfオブジェクトに対してst_transform関数を適用することで座標の投影変換を行うことができる。ここでは調査地点の位置を考慮し、「UTM Zone 52S」に投影変換を行っている。 投影返還を行った場合には座標単位が変わるため、coord_sfの範囲指定はそれに合わせた値を与える必要がある。ただし、描画時の軸ラベルは緯度経度表示のままとなる。
basemap_utm <- basemap %>%
st_transform(crs = "+proj=utm +zone=52 +south +datum=WGS84 +units=m +no_defs")
island_diversity_sf_utm <- island_diversity_sf %>%
st_transform(crs = "+proj=utm +zone=52 +south +datum=WGS84 +units=m +no_defs")
island_diversity_sf_utm
## Simple feature collection with 60 features and 15 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 673348.8 ymin: 9942260 xmax: 677390.6 ymax: 9944082
## CRS: +proj=utm +zone=52 +south +datum=WGS84 +units=m +no_defs
## # A tibble: 60 x 16
## island_ID n_species data data2 Shannon_H island_area island_perimeter
## * <chr> <int> <lis> <lis> <dbl> <dbl> <dbl>
## 1 GB1 18 <tib~ <tib~ 2.40 4774. 303.
## 2 GB10 5 <tib~ <tib~ 1.39 121. 46.4
## 3 GB11 9 <tib~ <tib~ 1.54 817. 120.
## 4 GB12 14 <tib~ <tib~ 2.08 1650. 154.
## 5 GB13 10 <tib~ <tib~ 2.12 602. 90.0
## 6 GB14 11 <tib~ <tib~ 1.99 535. 90.2
## 7 GB15 9 <tib~ <tib~ 1.58 381. 77.8
## 8 GB16 6 <tib~ <tib~ 1.20 137. 43.5
## 9 GB17 3 <tib~ <tib~ 0.509 18.4 18.1
## 10 GB18 8 <tib~ <tib~ 1.78 433. 77
## # ... with 50 more rows, and 9 more variables: distance_Gam <dbl>,
## # buffer_area_1000m <dbl>, tree_basal_area <dbl>, species_number <dbl>,
## # soil_depth_mean <dbl>, leaf_litter_cover <dbl>, no_transects <dbl>,
## # no_plots <dbl>, geometry <POINT [m]>
g_gis_utm <- ggplot() +
geom_sf(data = basemap_utm, color = "black", fill = "darkolivegreen3") +
geom_sf(data = island_diversity_sf_utm, aes(color = Shannon_H), size = 4, alpha = 0.8) +
scale_color_viridis_c(name = "Shannon's H", breaks = seq(0, 3, by = 1)) +
coord_sf(xlim = c(673000, 678000), ylim = c(9942000, 9945000)) +
labs(caption = "Land data from OSM(https://osmdata.openstreetmap.de/data/land-polygons.html)") +
theme_example +
theme(
legend.position = c(0.15, 0.20)
)
ggsave(g_gis_utm, file = "g_gis_utm.png", width = 14, height = 8, dpi = 200)
g_gis_utm
道具は道具に過ぎないから、目的に則したものを選ぼう
もしも困ったことがあったら、
お役立ちページ
連絡先:m.eco.3793(アットマーク)gmail.com